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Present and planned dark matter detection experiments search for WIMP-induced nuclear recoils 
in poorly known background conditions. In this environment, the maximum gap statistical method 
provides a way of setting more sensitive cross section upper limits by incorporating known signal 
information. We give a recipe for the numerical calculation of upper limits for planned directional 
dark matter detection experiments, that will measure both recoil energy and angle, based on the 
gaps between events in two-dimensional phase space. 
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I. INTRODUCTION 



Dark matter comprises approximately 25% of the mass 
of the universe [l[ . The generic dark matter candidate is 
a weakly interacting massive particle (WIMP). If WIMPs 
are supersymmctric particles, the predicted mass is in the 
range of 10 to 10 4 GeV/c 2 , and the expected cross section 
lies in the range of 10~ 42 to 10~ 48 cm 2 Q. Many experi- 
ments seek to detect dark matter particles via their elas- 
tic scattering interactions with detector nuclei Re- 
cent measurements 0, H| limit the cross section to be 
less than approximately 5 x 10~ 44 cm 2 . Given the small 
size of the expected WIMP cross section, discrimination 
against backgrounds is of paramount importance in di- 
rect detection experiments. For the same reason, there 
is much to gain by optimizing the statistical methods 
used to interpret experimental data as upper limits on 
the WIMP interaction cross section Q . 



In this paper, we develop a new statistic, the maximum 
patch, for setting limits on the WIMP-nucleus interac- 
tion cross section. This method is motivated by direc- 
tional dark matter experiments, which seek to measure 
both the nuclear recoil energy, and the recoil angle of the 
struck nucleon in WIMP-nucleon interactions. In section 
[II] we introduce the theoretical distributions used for set- 
ting limits, and in sections [TTT] and [TV] we discuss limit 
setting techniques within the context of one- and two- 
dimensional WIMP detection experiments. In section fVl 
we compare results in the cases of (i) observing a signal 
with no background, (ii) observing some signal and some 
background, and (iii) observing only background. 



II. SETTING DARK MATTER CROSS 
SECTION LIMITS 

Direct detection experiments typically measure the en- 
ergy deposited by the recoil nucleus [3|, infer the true 
nuclear recoil energy, and set upper limits on the WIMP- 
nucleus interaction cross section by comparing the the- 
oretical distribution with this one-dimensional data set. 
The theoretical event rate distribution is given by 0] 
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where E is the nuclear recoil energy, E = ^mr>v\ 
is the dark matter particle's kinetic energy, r = 
AmofTiT I '(mo + ?tjt) 2 with dark matter particle mass 
m D and target nucleus mass m T , v t hreshoid and v max 
are the minimum observable and escape velocities of the 
dark matter (taken to be the dark matter velocity that 
produces a maximum recoil energy equal to the exper- 
imental lower limit recoil energy threshold and oo here 
respectively, for simplicity), ve = 244 km/s is the Earth's 
velocity relative to the dark matter halo, and f(v, ve) is 
the dark matter velocity distribution function, assumed 
here to be a Maxwell-Boltzmann distribution with RMS 
velocity vo = 230 km/s. The normalization factor Rq 
is the event rate per unit mass with (vthreshoidiV max ) = 
(0,oo), 
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where Nq is Avogadro's number, A is the atomic mass of 
the target, pu is the dark matter density, taken here to be 
0.3 (GeV I c 2 ) / cm 3 , and <Jo is the zero momentum-transfer 
dark matter-nucleus interaction cross section. We use 
the fact that the differential interaction rate scales sim- 
ply with (To in the following discussion of limit setting 
techniques. 

A new thrust in the field of WIMP searches has been 
to develop detector sensitivity in a second dimension, the 
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nuclear recoil angle 

USES- The WIMP-nucleus inter- 
action signal is expected to be highly anisotropic in recoil 
direction because of the earth's motion with respect to 
the WIMP halo In contrast, the backgrounds of 

most WIMP experiments are relatively isotropic in recoil 
angle in the detector coordinate system [12], and there- 
fore this experimental approach can provide increased 
discrimination against backgrounds. It has recently been 
suggested that WIMP direct detection searches relying 
on a recoil energy signature alone may be insufficient for 
distinguishing WIMP events from nuclear recoils if the 
WIMP-nucleus cross section is smaller than the coher- 
ent scattering cross section for solar neutrinos [13$ . This 
makes directional detection particularly attractive, since 
solar neutrino-induced recoils point back to the sun, un- 
like recoils from WIMP interactions. Directional detec- 
tion could also potentially probe the velocity distribution 
of our galaxy's dark matter halo [l4|. The theoretical 
distribution as a function of nuclear recoil energy E and 
recoil angle ip (where ip is the angle in the detector lab 
frame between the nuclear recoil track observed in the de- 
tector and the direction the dark matter "wind" is blow- 
ing, which is normally taken to be the vector poin ting 
from the constellation Cygnus to Earth) is given by 

d 2 N(v E , oo) 1 Rp ( (vecosi/j - v min ) 2 \ 
dEd(cosip) ~ 2E r eXP \ v% J U 

where iVnm = (E ' j ' Eqt)^I 2 vq is the smallest dark matter 
particle velocity which can produce a nuclear recoil with 
energy E. 

Given a theoretical distribution, one can compare with 
an observation to set a limit on the WIMP-nucleus in- 
teraction cross section. The usual method to obtain an 
upper limit at some confidence level is to vary the theo- 
retical parameters until the appropriate cumulative prob- 
ability distribution function (CDF) takes on the confi- 
dence level desired (0.9 for a 90% confidence level up- 
per limit) when evaluated at the observed statistic (e.g. 
the number of observed events). In the following two 
sections, we discuss upper limit calculations using the 
one- and two-dimensional theoretical distributions re- 
spectively, together with several statistics of interest. 



III. DARK MATTER STATISTICS IN ONE 
DIMENSION 

Here we compare the traditional Poisson method with 
the maximum gap, a statistic often used in dark mat- 
ter experiments for obtaining an upper limit with one- 
dimensional data [1, Q. While the traditional Poisson 
method is based solely on the number of counts ob- 
served [l5j], the maximum gap procedure incorporates 
what is known about the shape of the expected signal 
into the limit determination [6|. For the discussion that 
follows, consider a series of nuclear recoil energy mea- 
surements {Ei, ...,En} where N is the total number of 



measurements. Assume the data points are distributed 
with a known theoretical function dN(X)/dE, where the 
A are the parameters of the theoretical model dN/dE 
given in equation [T] Assuming standard values for the 
dark matter halo parameters, there is only one free pa- 
rameter which we then vary to set upper limits, the 
zero-momentum transfer dark matter-nucleus interaction 
cross section ctq in equation^ 



A. Poisson Method 

A straightforward cross section upper limit on a given 
data set can be obtained by employing the Poisson 
method. To set a limit, we are interested in the prob- 
ability, given a value of the cross section <ro in a theoret- 
ical distribution dN/dE, that the total number of events 
observed in our data is equal to a certain value or less. 
If we are conservative and assume no knowledge of the 
background distribution and therefore that all observed 
events arc signal, then an upper limit at some desired 
confidence level may be set by adjusting cto in dN/dE 
until the total number of events \x expected, given by in- 
tegrating dN/dE over the whole experimental range, is 
such that it satisfies equation 0J 

a = E ^7 ( 4 ) 

Here, 1 — a is the confidence level of the upper limit set in 
this way, and N is the number of observed data events. In 
order to incorporate knowledge of expected backgrounds 
into equation 21 we must use a modified form of this rela- 
tion that assumes the overall normalization of the back- 
ground, which is often poorly understood in dark mat- 
ter direct detection experiments (for instance, see equa- 
tion 32.35 in [HI). The most conservative approach is 
to assume no knowledge of the backgrounds, and so all 
events arc signal candidates. In this case, any observed 
events considerably degrade the upper limit obtained 
with equation 3J This is particularly true for scenar- 
ios in which a background fills a small subset of the full 
experimental acceptance. For nuclear recoil signals, as 
the energy detection threshold is lowered, the sensitivity 
to backgrounds increases. Background events observed 
near detection threshold arc counted with the same sig- 
nificance as events in higher recoil energy sub-intervals of 
the experimental acceptance. Direct dark matter detec- 
tion experiments gain sensitivity to WIMP events by low- 
ering their energy thresholds, since the WIMP-nucleon 
interaction rate is expected to peak at low nuclear recoil 
energies. In this scenario, the Poisson method can lead 
to overly conservative upper limits on the WIMP-nucleon 
interaction cross section, in the presence of backgrounds. 



3 



B. Maximum Gap Method 

For a given <To, the "gap" for a pair of data points is 
defined to be M 



" El+1 dN 

— {a )dE 



(5) 



where X{ is the value obtained by integrating dN/dE 
between the observed energy values [Ei,E i+ i] for i = 
0,..,N (see figure 1 in and E and E N+1 are the 
lower and upper recoil energy experimental thresholds. 
A set of ./V recoil energy measurements yields an (N — 1)- 
dimensional vector x of gaps. The "maximum gap" for a 
set of N recoil energy measurements is defined to be the 
largest member of the set of all gaps that can be com- 
puted from the data. This quantity depends on an inte- 
gral over the hypothesized theoretical distribution, and 
through that integral, on the WIMP interaction cross 
section ctq ■ The larger the maximum gap given a ag , the 
larger the discrepancy between the number of points ob- 
served in data and the number of points expected. There- 
fore, the maximum gap allows a powerful statistical test 
between the measured data and the normalization of the 
theoretical signal distribution. 

To set a limit, we are interested in the probability, 
given a value of the cross section cr in a theoretical dis- 
tribution dN/dE, that the maximum gap for a given set 
of data is equal to a certain value or less. This is de- 
scribed by the CDF of the maximum gap, and can be 
analytically calculated [f|, 
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where p is the total number of events expected in the 
experimental range {dN / dE integrated from lower to up- 
per limit energy threshold), and m is the greatest integer 
< p/x. An upper limit at a given confidence level is 
obtained from equation [6] by adjusting the input cross 
section <to until the function Co above, evaluated at the 
maximum gap of the data, yields 0.90. The interpre- 
tation is that in an ensemble of experiments, on each of 
which the upper limit setting technique is employed, 90% 
will obtain an upper limit greater than <7o- 



C. Discussion of Limit Setting Techniques 

The maximum gap statistic possesses a number of 
nice qualities, particularly in the presence of background, 
which motivate the generalization of the method to two 
dimensions for directional dark matter detection experi- 
ments. We summarize the main results here; see |6| for 
a rigorous discussion. 

First, the maximum gap is unchanged under a one-one 
transformation of the variable in which the events are dis- 
tributed. This can be seen by making a transformation 



from recoil energy in the above discussion to a variable 
p such that p is equal to the total number of events ex- 
pected between the point E and the lower energy thresh- 
old. In p, equation[Tjis a uniformly distributed, unit den- 
sity function. This calculation is treated in more detail 
in Appendix I. This means that the maximum gap does 
not depend on the form of the theoretical distribution. 
Second, the method can be used for an arbitrarily large 
number of observed data points, and requires no binning 
of the data points. Most importantly, this statistic pro- 
vides a conservative upper limit on the true WIMP cross 
section that can considerably out-perform the Poisson 
upper limit setting technique. 

In WIMP detection experiments there arc often low 
energy backgrounds from processes which are difficult to 
model accurately with simulations, such as MeV-scale 
neutron interactions or 238 U and 232 Th decay progeny 
in the detector materials. It is unlikely that the mea- 
sured maximum gap will be found in an interval of 
an experiment's acceptance that contains both signal 
WIMP events and a large number of background events. 
The presence of sizable background would significantly 
shorten the gap sizes expected from signal alone. We 
thus expect the maximum gap to occur in the data in 
an interval where the background does not dominate. In 
this way, the maximum gap method automatically se- 
lects recoil energy intervals that are characteristic of the 
expected WIMP signal alone. 

Another striking advantage of the maximum gap 
method is that, for a fixed gap size, it is independent 
of the total number of events observed. In the Poisson 
case, each additional point observed inflates the upper 
limit an experimenter sets on their data. On the other 
hand, if an experimenter observes a large gap in their 
data, the limit set on that data is unchanged if the num- 
ber of points observed outside of the maximum gap is 
one or one million. 



IV. A NEW DARK MATTER STATISTIC FOR 
TWO DIMENSIONS 

For directional dark matter detection experiments, 
it is desirable to preserve the benefits of the maxi- 
mum gap method for setting upper limits. Towards 
this end, we need to generalize the method to two 
dimensions. We consider a scries of measurements 
cos(-0)i), {Em, cos(^)at)}, where N is the total 
number of measurements of the energy of nuclear recoils 
E, and ip is the nuclear recoil angle in a dark matter 
detection experiment measured from the vector pointing 
from the constellation Cygnus to the Earth in the detec- 
tor lab frame. In general, the two-dimensional rate will 
be given by a function d 2 N(X)/d(cos(tp))dE, where the 

A are the theoretical parameters of the model. As in one 
dimension, the model for the two-dimensional differential 
WIMP-nucleon interaction rate, equation[31 under stan- 
dard dark matter halo assumptions, depends only on ctq, 
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the true WIMP interaction cross section. In the follow- 
ing we describe an algorithm for obtaining general, multi- 
dimensional CDFs, focusing on the two-dimensional case. 
We then apply the usual prescription for setting upper 
limits, varying (To, and find that our two-dimensional 
limit setting technique has the correct 90% coverage. 



A. Monte Carlo Generated Cumulative 
Distribution Functions 

To calculate the CDFs for an arbitrary statistic of in- 
terest (SI), we simplify matters by asking an easier ques- 
tion than "what is the probability that the SI is less than 
or equal to a certain value" and instead ask "what is the 
probability that, given an observation of N events, the 
SI is less than or equal to a certain value." The advan- 
tage of the latter question is that it can be addressed 
with Monte Carlo methods, and leads naturally to the 
resolution of the first question. 

We start by generalizing the one-dimensional case. The 
theoretical distribution in equation [I] assuming a value 
for do, gives a concrete form for dN/dE, the expectation 
for how the observed events are distributed in nuclear 
recoil energy. Then, we draw N events from this dis- 
tribution. We compute the value of the SI on this fake 
data, whether it be the maximum gap of the distribution, 
or some other SI. We repeat this procedure many times, 
until we have an SI frequency distribution, given N ob- 
served events. We do this in turn for N = 1,2, N max 
stopping for N max so large that the Poisson probability 
(equation[H]) for observing N max events is negligible. This 
results in an array h — {ho, hN max } of histograms, 
where ho is the frequency plot of the SI given the obser- 
vation of zero events, h\ is the frequency plot of the SI 
given the observation of one event, etc. 

In each of the hi, i = 1, ...,N rnax , we have N t0 y S en- 
tries, where N t0 y S is the number of toy experiments, for 
each number of observed events, that we performed. By 
normalizing each of the hi by N t0 y S , we obtain the prob- 
ability distribution of the SI, for each of the i events 
observed subclasses considered. The resulting normal- 
ized vector of histograms can be properly interpreted as 
the probability distribution functions (PDFs) of the SI 
h = h/N toys . 

For setting an upper limit, we need the CDFs of the SI. 
To construct these, we take each histogram member hi of 
h in turn and create a new histogram, h c j, assigning to 
each bin b' rni of h c j the value given in equation where 
bki is the k th bin of hi, and k and m index the bins of 
the PDF and CDF histograms, respectively. 

Ill: 

fc=i 

For convenience, 500 bin CDF and PDF histograms were 
generated for the maximum gap studies in this paper, 



and 300 bin CDF and PDF histograms were generated 
for the maximum patch studies. 

In this way we numerically turn each of the PDFs in 
h into a CDF in h c for the SI. Now we have h c , a vector 
of CDFs, each corresponding to a number of observed 
events. The probability of observing a certain number 
of events is determined by the Poisson probability of ob- 
serving N events given jj, total events, 

u N 

P(M,N) = jfie-». (8) 

In order to construct the full CDF of the SI for the in- 
put theoretical distribution dN(ao)/dE, we add all of the 
histograms in h c together, weighting each by the prob- 
ability in equation [5] for seeing that number of events. 
This yields the following equation for the Monte Carlo 
generated CDF for the SI, 

C SI (x SI ,fx) = £ hkixsi) ( jje~» ) (9) 
fc=0 ^ ' ' 

Here xsi is the value of the SI at which we want to 
know the value of the CDF, and fi is the total number 
of events expected in the experimental range. The no- 
tation h Ct k(xsi) means to evaluate the value of the k th 
histogram in the h c vector of SI CDFs at xsi- This can 
be done in a number of different ways; we can take this as 
the height of the bin in h c ^ that xsi falls into (thereby 
obtaining a discrete CDF), or, perform a bin-to-bin in- 
terpolation, thus turning the h c ,k^ into smooth functions 
of xsi- In the results presented in this paper, the h c ^ 
histograms are interpolated using splines. 

We have now in equation [5] manufactured the analogue 
to equation [6] for the maximum gap statistic for an arbi- 
trary SI. Note however that for the maximum gap statis- 
tic this discussion is unnecessarily complicated because 
the maximum gap is unchanged under a one-one trans- 
formation of the theoretical distribution dN/dE in E, as 
shown in Appendix I. This property allows one to trans- 
form any distribution into a unit density, uniformly dis- 
tributed function. Thus whenever we change a® for the 
maximum gap statistic, we need not draw events from a 
different distribution, we may instead always use a unit 
density, uniform distribution. 

We validate the Monte Carlo CDF generating scheme 
by comparing the frequency distribution of upper limits 
resulting from the Monte Carlo version of the maximum 
gap method with the result obtained using the analytic 
CDF in equation [6j We find that the results agree within 
numerical errors. 



B. Maximum Patch Method 

Analogously to the one-dimensional gap, we may con- 
struct a "patch" as a subset of an experiment's total ac- 
ceptance, which is determined by (-Eqj cos(V>)o)j the en- 
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crgy and angle lower limit experimental threshold of mea- 
surement, and (-Ejv+i) C0S (' t P)N+i)i the energy and angle 
upper limit experimental threshold of measurement. Wc 
define a patch for a set of N data points to be 

Vijk = / -77 FmTp (fo)d£d(cos(V>)) 

JcoB(i>) h Je % d(cos{ip))dE 

(10) 

where i ranges from to N and j and k range, indepen- 
dently, from 1 to N. We require that cos(V')j > cos(V>)/c 
and that < £y < £? i+ i and J5j < E k < E i+1 . Wc 
also include the additional patch candidate not picked up 
by this prescription for every i; namely that one which 
has as its borders in cos(?/>) the upper and lower angular 
limits [cos(^>)o,cos(VO./v+i]- Equation [TU] describes rect- 
angular patches, whose limits are [£".;, E^+i] in E and 
[cos(V>)k,cos(^)j] (plus [cos(VOo,cos("0)jv+i]) in cos(V>). 
Some of the rectangles described by equation[Tn]may have 
points inside their boundaries; these are disqualified from 
being maximum patches, for the same reason that gaps 
with points in them are disqualified in the maximum gap 
method. In principle any two-dimensional shape can be 
used to define a patch; we have chosen to use rectangles 
for ease of computation. It is possible that sensitivity 
could be gained in this method by considering patch ge- 
ometries other than rectangles. 

To set an upper limit we are interested in the maxi- 
mum value that y^fc takes in equation [10] for all accept- 
able values of j. k and i. The situation is illustrated in 
figurc[l] where the spikes are Monte Carlo generated fake 
data points in E and cos(V>) (4 in total) and the smooth 
curve represents d 2 N (<jq) / d(cos(ip))dE for a given cross 
section value. The maximum patch candidate in this ex- 
ample has boundaries that extend from the lower to the 
upper cos(V0 threshold, and from 5—15 keV. Note that 
this particular maximum patch candidate is also a max- 
imum gap candidate in i?-space. 

Our algorithm for computing the maximum patch on a 
set of observed two-dimensional data points is described 
in detail in Appendix II. An example of the patch-finding 
algorithm for 1 observed event is shown in figured] Once 
the maximum patch is found, one calculates the PDFs 
and sums them, appropriately weighted, to produce the 
CDFs as in section IIV Al Figure [3] shows the resulting 
CDF for several expected numbers of events [i. 

We note that we are able to use an analogue to the 
simplified CDF generation scheme mentioned at the end 
of IIV Al for the maximum gap method, with one impor- 
tant change. Unlike the maximum gap statistic, the max- 
imum patch is not unchanged under a one-one transfor- 
mation of the theoretical distribution, d? N j d(cos(ip))dE , 
in E and cos(-0). Therefore, to set limits on data dis- 
tributed according to d 2 N/d(cos(ip))dE with CDFs gen- 
erated from a unit density, uniformly distributed func- 
tion, one must use the transformation given in Appendix 
C of 0. 




FIG. 1: An illustration of a maximum patch candidate. The 
spikes are hypothetical measured data points in an experi- 
ment, and the smooth curve is, for some assumed cross sec- 
tion, the expected distribution of signal between the two low- 
est energy data points. The volume under the curve is one 
of the yijk's of equation [10] for this dataset, and thus a maxi- 
mum patch candidate. Note that if we project the data points 
and theoretical distribution onto the energy interval, this is 
also a maximum gap candidate in the nuclear recoil energy 
dimension. 




FIG. 2: An example illustrating the maximum patch algo- 
rithm for 1 assumed data point in the experimental accep- 
tance. If the data point is not on one of the four bound- 
aries, there are four patches that contribute. Our algorithm 
is detailed in Appendix II; the order in which our algorithm 
computes the above maximum patch candidates is upper left, 
upper right, lower left and lower right. 



C. The Recipe 

The recipe for the experimenter wishing to use the 
maximum patch method to set an upper limit on the dark 
matter cross section in her experiment is as follows, as- 
suming a measurement of a vector of N two-dimensional 
data points D = {(Ei, cos(-0)i), (-Ejv> cos(^))jv)}. We 
take as an explicit example a direction sensitive dark 
matter direct detection experiment in this recipe, but this 
method can be used for any two-dimensional dataset for 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(patch size/|J,) 

FIG. 3: Cumulative probability distribution functions (CDFs) 
for various total expected numbers of events fi in the exper- 
imental acceptance (equation [9] for the maximum patch SI, 
equation I10|) . A horizontal line is drawn at the 90% confi- 
dence level. So, for instance, for /i = 5, 90% of the time, the 
maximum patch of a toy signal experiment will be less than 
w 4. 



which the distribution of the signal is known. 



1. Given a predicted two-dimensional rate 
d 2 N(ao)/d(cos(ilj))dE, the experimenter cal- 
culates the maximum patch of their data for some 
starting value <7i for the WIMP-nuclcon interaction 
cross section by applying the recipe for calculating 
the maximum patch of a set of two-dimensional 
data points outlined in Appendix II. 



D. Validation of the Maximum Patch Method 



We validate the maximum patch prescription described 
above by checking that the coverage of our upper limit 
setting technique is correct. To do this, we generate many 
ensembles of toy Monte Carlo experiments, each with dif- 
ferent WIMP-nucleon interaction cross sections, and with 
signal events distributed according to equation[3] We set 
an upper limit on each toy experiment using the max- 
imum patch technique. Note that to test coverage at 
the 90% confidence level, we must choose cross sections 
that yield an expected number of events /i > 2.3026. Be- 
low this, no cross section upper limit can be set with a 
confidence level as high as 90% (see figure [3]). For each 
input cross section <ro, and the associated total number of 
events /i, we perform 10,000 toy experiments, where the 
observed number of events is Poisson-distributcd about 
/U, and (E,cos(i/j)) for the events are distributed accord- 
ing to equation [3] For each ensemble of 10,000 experi- 
ments with different input cto's, we record the percentage 
of the time our upper limit is higher than co- For an up- 
per limit requested at 90% CL, the percentage should 
be 90% within statistical errors, which is the definition 
of correct coverage. The results of this study, for era's 
such that /i = 1, ...,30 are shown in figure [4] The step- 
ping behaviour observed in the Poisson limit coverage 
is expected, due to the discrete nature of the statistic. 
Figure 0] also shows the coverage computed in this way 
for the maximum gap method (with the CDFs computed 
via our Monte Carlo technique) and the Poisson method. 
All methods are observed to have the correct coverage, 
within statistical and numerical errors. 



2. The experimenter then must evaluate equation [5] 
for the case where the statistic of interest is the 
maximum patch. The maximum patch CDF for /i 
total expected events can either be calculated by 
following the Monte Carlo procedure outlined in 
section HV Al or by referencing the tables provided 
at the end of this paper in Appendix III. 

3. Evaluating equation[9]at the maximum patch of the 
data yields the confidence level at which u\ is an 
upper limit on the WIMP-nucleon interaction cross 
section. If this is less (more) than the confidence 
level desired, the cross section guess is increased 
(decreased) to a new guess <T2. The experimenter 
again calculates the maximum patch of the data for 
this assumed cross section and then, via equation[9j 
calculates the CL at which <T2 is an upper limit on 
the WIMP-nucleon interaction cross section. This 
procedure is iterated for as many guesses <jn as 
required to set the desired confidence level upper 
limit on the WIMP-nucleon cross section for the 
data. 
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FIG. 4: For toy Monte Carlo experiments with pure signal 
events, the coverage as a function of signal events for the 
maximum gap, maximum patch and Poisson methods. This 
demonstrates that our upper limit setting methods have the 
correct coverage, within statistical errors for /i > 2.23026, as 
expected. Each point in this plot corresponds to 10,000 MC 
toy experiments. The errors shown are statistical. 
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V. COMPARISON OF METHODS 

Having built the maximum patch method, we compare 
it to various other methods for setting upper limits, in 
several circumstances of interest. Our goal is to highlight 
the impact of directionality for dark matter direct detec- 
tion experiments. Unless otherwise stated, for the various 
comparisons in this paper we arbitrarily assume a WIMP 
mass of 60 GeV, and we use the XenonlO experiment's 
acceptance and target properties [3] to construct limits 
(ignoring subtleties like quenching factors, form factors 
and efficiencies). 

First, in the absence of background and with a sizable 
signal, we would like to verify that the maximum patch 
method has not only the same coverage as the Poisson 
method (see subsection HVDp but also that its perfor- 
mance is comparable as a method for setting upper limits. 
Towards this end, we employ the ensembles of 10,000 toy 
Monte Carlo experiments from subsection IIVD1 record- 
ing the median upper limit obtained by the maximum 
patch, maximum gap and Poisson methods as a function 
of /i for each 10,000 event sub-sample generated with 
a different input a^. This comparison is shown in fig- 
ure O Figure [B] shows the frequency distribution of up- 
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FIG. 6: The frequency distribution of upper limits auL ob- 
tained by applying the Poisson, maximum gap and maximum 
patch procedures to 10,000 toy Monte Carlo experiments with 
an input cross section <7o shown as a black line on the graph. 
o"o was chosen to yield 7 total expected events in the exper- 
imental interval. The percent of the time that our proce- 
dure for each method sets a correct upper limit (a limit above 
the true input cross section oo) is (90 ± 1)%, (90 ± 1)% and 
(92 ± 1)% for the maximum gap, maximum patch and Poisson 
methods, respectively. 



-e- poisson 
-B max gap 
■ a max patch 




5 10 15 20 25 30 

number of signal events 

FIG. 5: The median upper limit cross section obtained by our 
implementation of the Poisson, maximum gap and maximum 
patch procedures divided by the true input cross section exo as 
a function of input cross section (and thereby, as a function 
of the total number of expected events in the experimental 
interval). Each point in this plot corresponds to 10,000 MC 
toy experiments. 

per limits set by the maximum gap, maximum patch and 
Poisson methods on the 10,000 event ensemble with an 
input do such that a total of /i = 7 signal events are ex- 
pected. These three distributions are used to generate 
the fx = 7 point in both figure [5] and figure 2J The com- 
puted coverages in figure [6] for the maximum gap, maxi- 
mum patch and Poisson methods are (90±1)%, (90±1)% 
and (92 ± 1)%, respectively, which is correct, within sta- 
tistical errors. 

We note that in the case of pure signal, the Poisson 
method outperforms the maximum gap and maximum 



patch methods by a factor of ~ 1.2. This is expected; the 
maximum gap procedure has already been demonstrated 
to give looser upper limits than the Poisson method in 
the case of pure signal (see @, figure 3(a)). This could 
be resolved, as in [|| by considering gaps containing 
greater than zero events, which is termed the optimal 
gap method. The extension of the optimal gap approach 
to two dimensions is not considered here. We also note 
that little sensitivity is gained by using the maximum 
patch method versus the maximum gap method in the 
case of pure signal. 

Realistically, WIMP direct detection experiments have 
backgrounds, and so a more interesting comparison is 
how the maximum gap, maximum patch and Poisson 
methods do in the presence of sizable backgrounds. For 
this test, we populate the lower half of our nuclear recoil 
energy range, and the lower half of our nuclear recoil an- 
gular range, with a background drawn according to a flat 
distribution. We caution that all of our results including 
background events are highly dependent on this partic- 
ular background distribution choice. Figure [7] shows the 
frequency distribution of upper limits obtained by ap- 
plying the maximum gap method, the maximum patch 
method, and the Poisson method to 10,000 toy exper- 
iments generated in this way with a total number of 
expected background events of 7, and a total expected 
number of WIMP signal events of 5. The coverage is 
(100 ± 1)%, (95 ± 1)% and (100 ± 1)% for the maximum 
gap, maximum patch and Poisson methods respectively, 
which is not at the confidence level requested due to the 
presence of the large background. The median 90% con- 
fidence level upper limit cross sections, from the maxi- 
mum gap, maximum patch and Poisson techniques are 
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compared in figure H as a function of the total num- 
ber of expected input background events [i = 1, 2, 30 
(with 5 expected signal events in each toy experiment). 
The total number of background events in a given Monte 
Carlo experiment, like the total number of signal events, 
is drawn randomly from a Poisson distribution with the 
mean given by the total number of expected events. 

Figure [8] shows that in the presence of a sizable WIMP 
signal, the maximum patch procedure provides stronger 
upper limits than the Poisson or maximum gap proce- 
dures as the amount of background contamination in- 
creases. The Poisson method does so poorly because it 
uses only the total number of events to set upper limits, 
assuming that they are all signal, yielding an increasingly 
inflated upper limit as the number of background events 
injected into the toy experiments increases. The maxi- 
mum patch method outperforms the maximum gap pro- 
cedure because it includes an extra dimension in which 
signal and background are differently distributed. The 
observation that the maximum gap and maximum patch 
limits seem to asymptotically flatten is due to the overly 
simplistic background chosen for these studies; eventu- 
ally the maximum patch or gap is always outside of the 
lower E interval or E — cos(ip) quadrant, and thus char- 
acteristic only of the WIMP-signal input which, within 
statistical fluctuations, is identical outside of the lower E 
interval or E — cos(V>) quadrant as the backgrounds are 
confined there. 
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FIG. 7: The frequency distribution of upper limits (Jul ob- 
tained by applying the Poisson, maximum gap and maximum 
patch procedures to 10,000 toy Monte Carlo experiments with 
an input cross section <to shown as a black line on the graph. 
o"o was chosen to yield 5 total expected signal events in the 
experimental interval. Background events were included in 
these experiments with a uniform distribution in the lower 
half of the recoil energy and recoil angle experimental accep- 
tance. For this plot, the total number of expected background 
events was set at 7. The percent of the time that our proce- 
dure for each method sets a correct upper limit (a limit above 
the true input cross section oo) is (100 ± 1)%, (95 ± 1)% and 
(100 ± 1)% for the maximum gap, maximum patch and Pois- 
son methods, respectively. 

The most probable situation in current dark matter 




5 10 15 20 25 30 

number of background events 



FIG. 8: The median upper limit cross section obtained by our 
implementation of the Poisson, maximum gap and maximum 
patch procedures divided by the true input cross section ao as 
a function of background events injected into our toy exper- 
iments, uniformly distributed in the lower half recoil energy 
and recoil angle interval. The point at number of background 
events equals 7 comes from dividing the medians of the fre- 
quency distributions in figure[7]by the true input cross section 
for the Poisson, maximum gap and maximum patch limit set- 
ting procedures. Each point in this plot corresponds to 10,000 
MC toy experiments and has, in addition to background, a sig- 
nal component generated with a cross section corresponding 
to 5 total expected signal events. 



experiments is that the true cross section lies well below 
the detectable range. To study this scenario, we perform 
10,000 toy experiments as above, with oo = 1 x 10~ 46 
cm 2 (or « total expected signal events), and the same 
distributions of background events as in figures [7] and 
[8] (flat in the E — cos(-0) plane, and confined to the 
lower E — cos(-0) quadrant). Figure [9] shows the fre- 
quency distribution of upper limits obtained by applying 
the maximum gap method, the maximum patch method, 
and the Poisson method to 10,000 toy experiments with 
an expected number of background events of 7 and neg- 
ligible signal. The coverage is (100 ± 1)%, (100 ± 1)% 
and (100 ± 1)% for the maximum gap, maximum patch 
and Poisson methods, respectively, which is not at the 
confidence level requested due to the presence of a large 
background and no signal. The median 90% confidence 
level upper limit cross sections, from the maximum gap, 
maximum patch and Poisson techniques are compared in 
figure [Tol as a function of input background events in the 
case of negligible signal. 

In the presence of a negligible WIMP signal and in- 
creasing backgrounds, the maximum patch procedure 
provides by far the most restrictive upper limit as the 
amount of background contamination increases. From 
figure [10] it is clear that the Poisson method limit is not 
competitive as the number of backgrounds increases, and 
the maximum patch method outperforms the maximum 
gap method by at least a factor of 2 for more than 1 ex- 
pected background events for this particular background 
distribution. 
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FIG. 9: The frequency distribution of upper limits auL ob- 
tained by applying the Poisson, maximum gap and maximum 
patch procedures to 10,000 toy Monte Carlo experiments with 
an input cross section oo = 1 X 10~ 46 shown as a black line 
on the graph, chosen to give « events for the Monte Carlo 
toy exposures. Background events were included in these ex- 
periments with a uniform distribution in the lower half of the 
recoil energy and recoil angle experimental acceptance. For 
this plot, the total number of expected background events was 
set at 7. The percent of the time that our procedure for each 
method sets a correct upper limit (a limit above the true in- 
put cross section a ) is (100 ±1)%, (100 ±1)% and (100 ±1)% 
for the maximum gap, maximum patch and Poisson methods, 
respectively. 



In typical dark matter detection experiments, upper 
limits on the WIMP-nucleon cross section are reported 
as a function of WIMP mass. In order to compare the 
performance of the maximum gap, maximum patch and 
Poisson methods of setting upper limits as a function 
of WIMP mass, we generated 10,000 Monte Carlo toy 
datasets at a variety of different WIMP masses. Adopt- 
ing the likely scenario for a given dark matter experi- 
ment, we hypothesize a WIMP-nucleon interaction cross 
section of <jq = 1 X 10 -46 for these toy experiments and 
a background with an average number of events equal 
to 10 (again, flat in the E — cos(^) plane, and confined 
to the lower E — cos(i}i) quadrant). The result of this 
study is the three limit curves in figure QT] The solid 
bands represent the RMS widths of the upper limit fre- 
quency plots (much like figures [5J [7] and [5]) obtained on 
the 10,000 experiments at each given mass point. The 
shape of the maximum gap limits at low WIMP masses 
is an artifact of the background choice. For low WIMP 
masses, the maximum patch method is by far the most 
sensitive to the cross section (between w 10 — 100 GeV) 
with a very small RMS. At higher WIMP masses, the dif- 
ference between the maximum gap and maximum patch 
limit techniques diminishes, but both consistently out- 
perform the Poisson limit setting method. We note that 
the studies in this paper are all based on one assumed 
background shape. Different background distributions 
will lead to different results. 
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FIG. 10: The median upper limit cross section obtained by 
our implementation of the Poisson, maximum gap and max- 
imum patch procedures divided by the true input cross sec- 
tion (To as a function of background events injected into our 
toy experiments, uniformly distributed in the lower half recoil 
energy and recoil angle interval for a negligible input signal 
cross section. The point at number of background events 
equals 7 comes from dividing the medians of the frequency 
distributions in figure [5] by the true input cross section for 
the Poisson, maximum gap and maximum patch limit setting 
procedures. Each point in this plot corresponds to 10,000 MC 
toy experiments. 
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FIG. 11: The median upper limit cross section obtained by 
our implementation of the Poisson, maximum gap and maxi- 
mum patch procedures divided by the true input cross section 
o"o as a function of assumed WIMP mass, with 10 background 
events on average that are uniformly distributed in the lower 
half recoil energy and recoil angle interval for a negligible in- 
put signal cross section. Each point in this plot corresponds 
to 10,000 MC toy experiments. The bands correspond to the 
RMS widths of the upper limit frequency distributions ob- 
tained using the Poisson, maximum gap and maximum patch 
methods, respectively, at each mass point sampled. 
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VI. CONCLUSIONS 

In this paper, we have developed a new method for 
setting upper limits in two dimensions. The motivation 
for our maximum patch method is directional dark mat- 
ter detection, but it is generally applicable to any two- 
dimensional data sets for which the distribution of the 
signal is known. The approach is an extension of the one- 
dimensional maximum gap method [f|, a statistic often 
used to set limits on the WIMP-nucleon interaction cross 
section in direct detection dark matter experiments. To 
directly detect dark matter requires unprecedented con- 
trol and understanding of backgrounds. The great ad- 
vantage of the maximum gap and patch methods is that 
they require no knowlege of the background distribution 
to set conservative upper limits on the WIMP nucleon 
scattering cross section. The scattering kinematics of a 
dark matter signal in one and two-dimensional direct de- 
tection experiments are relatively simple. This informa- 
tion is included in a straightforward way in the maxi- 
mum gap and patch methods. In particular, the maxi- 
mum patch method incorporates information about the 
large expected angular anisotropy of dark matter scat- 
tering into the limit setting procedure. We demonstrate 
that for simplistic background assumptions, the maxi- 
mum patch method and two-dimensional dark matter de- 
tection yield a large gain in sensitivity, especially at low 
WIMP masses, over one-dimensional dark matter detec- 
tion. 



Appendix I 

Our goal is to prove, in detail, the assertion in [(J that 
the maximum gap is unchanged under a one-one transfor- 
mation of the variable in which the events are distributed, 
and thus that the maximum gap is independent of the 
particular way in which the events arc distributed for a 
given (To m one dimension (by dN (a®) / dE) . The total 
number of events in the experimental window is given by 
p, where 



a,nd dE /dp— (dp/ dE) 1 . dp/dE from equation[T2lis then 
given by the fundamental theorem of calculus to be 



+ ™dE d JL. 

dE 



(11) 



Here, if E is outside of the experimental thresholds, 
dN/dE = 0, and dN/dE is strictly positive. 

We wish to change variables from E in equation QT] to 
p, where p @ is 



By the chain rule, 



dN _ dN dE 
dp dE dp 



(12) 



(13) 




dE dE \ 7_ 00 dE' I dE 
Thus, equation 1131 implies that 
dN dNfdpY 1 



dp dE \dE 



> dE , 



115) 



and that in the new variable p, the differential rate is a 
unit density function, of length p. 



Appendix II 

Our method for calculating the maximum patch is 
as follows. A similar scheme is put forth in p^ |. 
Call the set of measured two-dimensional data points 
on which we wish to determine the maximum patch, 
D = {(Ei, cos('0)i), {En, cos(tp) n)} , where as above, 
(E ,cos(ip)o) and (-Ejv+i, cos^jv+i) wc define to be the 
upper and lower recoil energy and recoil angle thresholds 
of our experiment. 

1. Order D in E. Call the new, ^-ordered vector D'. 

2. Loop through the intervals [Ei, Ej) in D' such that 
i runs from to N + 1 and j runs from i + 1 to 
N + l. 

3. Inside each E interval in this loop, if there are no 
points with an Ek such that Ei < Ek < Ej, then 
compute equation 1101 with E limits [Ei,Ej] and 
cos(V>) limits [cos(^)o, cos(i/>);v+i]- 



4. If there are N p points with an Ek such 
that Ei < Ek < Ej, make a list of 
them, P = {(Eo,cos(ip)o),...,{E k ,cos(ip)k),..., 
(En + i,cos(iP)n+i)} of length N p + 2, ordered 
in E, adding the points (Eo, cos('0)o) and 
(-E/v+i, cos(^)jv+i) (the threshold points - that's 
where the +2 comes from in the total number of 
points in P) to the front, and back of the list, 
respectively. Order P by E, calling the new E- 
ordcred vector P' . Then loop through the inter- 
vals [cos(-0)_P',m, cos(V')p',n], in P' such that m 
runs from to Np + 2 and n runs from m + 1 to 
N p + 2, where the subscript (P',m) denotes the 
m th member of the ordered vector P' . Each iter- 
ation of this loop will provide a maximum patch 
candidate with E limits [Ei,Ej] and cos(-0) limits 

[cOS(V>)p',m,COs(V>)F',n]- 

5. For each of these candidates, check to see if there 
is a point D' d in D' such that Ei < Ed < Ej and 
cos(V0p',to < cos(V')d < cos(?/>)p',n- If so, then 
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there is a point inside our patch candidate, which 
means it is not a maximum patch candidate; throw 
it out. 

6. If a maximum patch candidate has passed all of 
the above criterion, stick it in a vector of some ar- 
bitrary length y. This vector exhausts all possible 
rectangles in the two-dimensional E-cos(tjj) plane. 

7. To find the maximum patch of the data, loop 
through all of the elements yi of y, and store the 
largest yi value; this is the maximum patch. 

Appendix III 

Tables I and II below record the h C! k{xsi)'s of the max- 
imum patch statistic for k = 1, 100 {h C fi(xsi) is for 
all maximum patch values except for the patch equal to 
the total number of expected events). The maximum 
patch CDF is computed by interpolating these points into 
smooth curves and using them to evaluate equation [9] 
Note that the smoothed curves obtained from the tables 
are functions of the maximum patch value divided 

by the total number of expected events. The work in this 
paper was done using 200 CDFs in the sum of equation[9] 



In Tables I and II, a deviation from a cumulative proba- 
bility of 100% for (y/fi) = 1 is observed at the 0.1% level 
for (j, = 67 and greater. The reader is advised to use the 
tables below only for \i = 50 or less, where they have 
been verified to give correct coverage at the 1% level. 

The reader is cautioned that these tables were gen- 
erated using a unit density, uniformly distributed the- 
oretical function as input, and thus cannot be directly 
applied, as in the maximum gap case, to data to obtain 
an upper limit. The data must be transformed such that 
it is uniformly distributed assuming a given model as 
its true distribution. The necessary transformation can 
be found in Appendix C of [l6l |. This subtlety can be 
avoided by generating a set of CDFs for every param- 
eter change in the theoretical model considered for the 
data, but in most applications this will be prohibitively 
computationally intensive. 
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TABLE I: Table of maximum patch CDF's (the h c ,k(xsi)'s of section H V Afl given an observation of n — 1, 54 events. These 
values of the CDFs are given as a function of the observed maximum patch divided by the total expected number of events. 
Upper limits can be set on two-dimensional data with these tabulated CDFs following the recipe laid out in section HV CI 
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TABLE II: Table of maximum patch CDF's (the h Ct k(xsi)'s of section TlV A|) given an observation of n = 55, 100 events. 
These values of the CDFs are given as a function of the observed maximum patch divided by the total expected number of 
events. Upper limits can be set on two-dimensional data with these tabulated CDFs following the recipe laid out in sect ion HV CI 
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